
#include "testHere.h"

double trace2(int m, int n, double *A, int lda, double *B, int ldb){
	double trace=0;
	int i;
	for(i=0;i<m;i++)
		trace += ddot(n,&A[i],lda,&B[i*ldb],1);
	return trace;
}

double myTrace2(char transa,char transb, int m, int n, double *A, int lda, double *B, int ldb){
	double trace=0;
	int i;
	if(transa=='N' || transa== 'n'){
		if(transb=='N' || transb== 'n')
			for(i=0;i<m;i++)
				trace += ddot(n,&A[i],lda,&B[i*ldb],1);
		else if(transb=='T' || transb== 't')
			return ddot(m*n,A,1,B,1);
			//for(i=0;i<m;i++)
			//	trace += ddot(n,&A[i],lda,&B[i],ldb);
	}
	else if(transa=='T' || transa== 't'){
		if(transb=='N' || transb== 'n')
			return ddot(m*n,A,1,B,1);
			//for(i=0;i<m;i++)
			//	trace += ddot(n,&A[i*lda],1,&B[i*ldb],1);
		else if(transb=='T' || transb== 't')
			for(i=0;i<m;i++)
				trace += ddot(n,&A[i*lda],1,&B[i],ldb);
	}
	return trace;
}

void testMyTrace2(){
	double *A, *B;
	int m=6,n=10;
	A = (double *) malloc( sizeof(double) * m * n);
	B = (double *) malloc( sizeof(double) * n * m);
	acml_randu(A,m*n,0,1);
	acml_randu(B,m*n,0,1);
	printf("\nThe matrix A\n");
	printM(A,m,n);
	printf("\nThe matrix B\n");
	printM(B,n,m);

	printf("The trace of A*B is %lf\n",myTrace2('N','N',m,n,A,m,B,n));
	printf("The trace of B*A is %lf\n",myTrace2('N','N',n,m,B,n,A,m));
	printf("The trace of BT*AT is %lf\n",myTrace2('T','T',m,n,B,n,A,m));
	printf("The trace of AT*BT is %lf\n",myTrace2('T','T',n,m,A,m,B,n));

	acml_randu(A,m*n,0,1);
	acml_randu(B,m*n,0,1);
	printf("\nThe matrix A\n");
	printM(A,m,n);
	printf("\nThe matrix B\n");
	printM(B,m,n);
	printf("The trace of AT*B is %lf\n",myTrace2('T','N',n,m,A,m,B,m));
	printf("The trace of B*AT is %lf\n",myTrace2('N','T',m,n,B,m,A,m));
	printf("The trace of A*BT is %lf\n",myTrace2('N','T',m,n,A,m,B,m));
	printf("The trace of BT*A is %lf\n",myTrace2('T','N',n,m,B,m,A,m));

}

/* 
* A is m*n B is n*k, C is k*m so final square matrix will be m*m whose trace is returned
*/
double myTrace3(char transa, int m, int n, int k, double* A, int lda, double* B, int ldb, double* C, int ldc){
	double *tmp;
	double trace;
	tmp = (double *) malloc( sizeof(double) * n*m);
	dgemm('N','N',n,m,k,1,B,ldb,C,ldc,0,tmp,n);
	trace = myTrace2(transa,'N',m,n,A,lda,tmp,n);
	free(tmp);
	return trace;
}

void testMyTrace3(){
	double *A, *B, *C;
	int m=6,n=10, k=12;
	A = (double *) malloc( sizeof(double) * m * n);
	B = (double *) malloc( sizeof(double) * n * k);
	C = (double *) malloc( sizeof(double) * k * m);
	acml_randu(A,m*n,0,1);
	acml_randu(B,n*k,0,1);
	acml_randu(C,k*m,0,1);
	printf("\nThe matrix A\n");
	printM(A,m,n);
	printf("\nThe matrix B\n");
	printM(B,n,k);
	printf("\nThe matrix C\n");
	printM(C,k,m);

	printf("The trace of A*B*C is %lf\n",myTrace3('N',m,n,k,A,m,B,n,C,k));
	printf("The trace of C*A*B is %lf\n",myTrace3('N',k,m,n,C,k,A,m,B,n));
	printf("The trace of B*C*A is %lf\n",myTrace3('N',n,k,m,B,n,C,k,A,m));
}


double slepNuclearNorm(int m, int n, double *A, int lda){
	double *sing,  info, norm, k=minof2(m,n);
	sing = (double *)malloc(sizeof(double)*k );
	dgesvd('N','N',m,n,A,lda,sing,NULL,1,NULL,1,&info);
	for (i=0; i<k;i++)
		norm += sing[i];
	free(sing);
	return norm;
}